/*
 * Matrix.java
 *
 * Created on July 18, 2006, 3:39 PM
 *
 */

package galronnlp.math;

import java.lang.Math;
import java.util.Vector;

/**
 *
 * @author Daniel A. Galron
 */
public class Matrix {

    private int row;
    private int col;
    private Double det;
    private Vector<Vector<Double>> matrix;
    
    private Matrix LU = null;
    private int[] pivotArray;
    /** Creates a new instance of Matrix */
    public Matrix(int row, int col) {
        this.row = row;
        this.col = col;
        matrix = new Vector<Vector<Double>>();
        for(int i = 0; i < row; i++) {
            Vector<Double> r = new Vector<Double>();
            //r.setSize(col);
            for(int j = 0; j < col; j++) {
                r.add(new Double(0.0));
                //r.set(j, new Double(0.0));
            }
            matrix.add(r);
        }
    }
    public Matrix transpose() {
        Matrix ret = new Matrix(col, row);
        for(int i = 0; i < row; i++) {
            for(int j = 0; j < col; j++) {
                try {
                    ret.put(j, i, this.get(i, j));
                } catch (MatrixIndexException ex) {
                    ex.printStackTrace();
                }
            }
        }
        return ret;
    }
    public double get(int r, int c) throws MatrixIndexException {
        if(r < 0 || r >= row || c < 0 || c >= col)
            throw new MatrixIndexException(r, c);
        return matrix.get(r).get(c).doubleValue();
    }
    public void put(int r, int c, double v) throws MatrixIndexException {
        if(r < 0 || r >= row || c < 0 || c >= col)
            throw new MatrixIndexException(r, c);
        matrix.get(r).set(c, new Double(v));
    }
    public boolean equals(Object o) {
        Matrix other = (Matrix)o;
        if(this.row != other.row || this.col != other.col) {
            return false;
        } else {
            for(int i = 0; i < row; i++) {
                for(int j = 0; j < col; j++) {
                    try {
                        if(this.get(i, j) != other.get(i, j))
                            return false;
                    } catch (MatrixIndexException ex) {
                        ex.printStackTrace();
                    }
                }
            }
            return true;
        }
    }
    /**
     * Matrix multiplication
     */
    public Matrix mult(Matrix b) throws MatrixSizeException {
       if(this.col != b.row) 
            throw new MatrixSizeException(this.row, this.col, b.row, b.col);
        Matrix ret = new Matrix(this.row, b.col);
        for(int i = 0; i < this.row; i++) {
            for(int j = 0; j < b.col; j++) {
                for(int k = 0; k < this.col; k++) {
                    try {
                        ret.put(i, j, ret.get(i, j) + (this.get(i, k) * b.get(k, j)));
                    } catch (MatrixIndexException ex) {
                        ex.printStackTrace();
                    }
                }
            }
        }
        return ret;
    }
    /**
     * Scalar multiplication (double)
     */
    public Matrix mult(double s) {
        Matrix ret = new Matrix(this.row, this.col);
        for(int i = 0; i < this.row; i++) {
            for(int j = 0; j < this.col; j++) {
                try {
                    ret.put(i, j, s * this.get(i, j));
                } catch (MatrixIndexException ex) {
                    ex.printStackTrace();
                }
            }
        }
        return ret;
    }
    /**
     * Scalar multiplication (int)
     */
    public Matrix mult(int s) {
        Matrix ret = new Matrix(this.row, this.col);
        for(int i = 0; i < this.row; i++) {
            for(int j = 0; j < this.col; j++) {
                try {
                    ret.put(i, j, s * this.get(i, j));
                } catch (MatrixIndexException ex) {
                    ex.printStackTrace();
                }
            }
        }
        return ret;
    }
    /**
     * Scalar multiplication (float)
     */
    public Matrix mult(float s) {
        Matrix ret = new Matrix(this.row, this.col);
        for(int i = 0; i < this.row; i++) {
            for(int j = 0; j < this.col; j++) {
                try {
                    ret.put(i, j, s * this.get(i, j));
                } catch (MatrixIndexException ex) {
                    ex.printStackTrace();
                }
            }
        }
        return ret;
    }
    /**
     * Matrix addition
     */
    public Matrix add(Matrix b) throws MatrixSizeException {
        if(this.row != b.row || this.col != b.col)
            throw new MatrixSizeException(this.row, this.col, b.row, b.col);
        Matrix ret = new Matrix(this.row, this.col);
        for(int i = 0; i < this.row; i++) {
            for(int j = 0; j < this.col; j++) {
                try {
                    ret.put(i, j, this.get(i, j) + b.get(i, j));
                } catch (MatrixIndexException ex) {
                    ex.printStackTrace();
                }
            }
        }
        return ret;
    }
    /**
     * Matrix substraction
     */
    public Matrix sub(Matrix b) throws MatrixSizeException {
        if(this.row != b.row || this.col != b.col)
            throw new MatrixSizeException(this.row, this.col, b.row, b.col);
        Matrix ret = new Matrix(this.row, this.col);
        for(int i = 0; i < this.row; i++) {
            for(int j = 0; j < this.col; j++) {
                try {
                    ret.put(i, j, this.get(i, j) - b.get(i, j));
                } catch (MatrixIndexException ex) {
                    ex.printStackTrace();
                }
            }
        }
        return ret;
    }
    public Matrix copy() {
        Matrix ret = new Matrix(this.row, this.col);
        for(int i = 0; i < this.row; i++) {
            for(int j = 0; j < this.col; j++) {
                try {
                    ret.put(i, j, this.get(i, j));
                } catch (MatrixIndexException ex) {
                    ex.printStackTrace();
                }
            }
        }
        return ret;
    }
    public String toString() {
        String ret = "[";
        for(int i = 0; i < this.row; i++) {
            ret += "[";
            try {
                ret += this.get(i, 0);
            } catch (MatrixIndexException ex) {
                ex.printStackTrace();
            }
            for(int j = 1; j < this.col; j++) {
                try {
                    ret += ", " + this.get(i, j);
                } catch (MatrixIndexException ex) {
                    ex.printStackTrace();
                }
            }
            ret += "]\n";
        }
        ret += "]\n";
        return ret;
    }
    
    public Matrix getSubMatrix(int i0, int j0, int i1, int j1) {
        Matrix ret = new Matrix(i1-i0+1, j1-j0+1);
        for(int i = i0; i < i1+1; i++) {
            for(int j = j0; j < j1+1; j++) {
                try {
                    ret.put(i-i0, j-j0, this.get(i, j));
                } catch (MatrixIndexException ex) {
                    ex.printStackTrace();
                }
            }
        }
        return ret;
    }
    
    public Matrix getSubMatrix(int[] r, int[] c) {
        Matrix ret = new Matrix(r.length, c.length);
        for(int i = 0; i < r.length; i++) {
            for(int j = 0; j < c.length; j++) {
                try {
                    ret.put(i, j, this.get(r[i], c[j]));
                } catch (MatrixIndexException ex) {
                    ex.printStackTrace();
                }
            }
        }
        return ret;
    }
    
    public Matrix getSubMatrix(int i0, int i1, int[] c) {
        Matrix ret = new Matrix(i1-i0+1, c.length);
        for(int i = i0; i < i1+1; i++) {
            for(int j = 0; j < c.length; j++) {
                try {
                    ret.put(i-i0, j, this.get(i, c[j]));
                } catch (MatrixIndexException ex) {
                    ex.printStackTrace();
                }
            }
        }
        return ret;
    }
    
    public Matrix getSubMatrix(int[] r, int j0, int j1) {
        Matrix ret = new Matrix(r.length, j1-j0+1);
        for(int i = 0; i < r.length; i++) {
            for(int j = j0; j < j1+1; j++) {
                try {
                    ret.put(i, j-j0, this.get(r[i], j));
                } catch (MatrixIndexException ex) {
                    ex.printStackTrace();
                }
            }
        }
        return ret;
    }
    
    // This is a re-implementation of the LU Decomposition algorithm as implemented
    // in the JAMA package. I modified it to conform to the internal structure of
    // this Matrix object (i.e. Vectors of Vectors of Double, rather than double[][])
    private void LUDecomposition() {
        if(LU == null) {
            LU = this.copy();
            pivotArray = new int[row];
            for(int i = 0; i < row; i++) {
                pivotArray[i] = i;
            }
            int pivsign = 1;
        
            Vector<Double> LUrowi;
            Vector<Double> LUcolj = new Vector<Double>(row);
        
            // Outer loop
            for(int j = 0; j < col; j++) {
                // Make a copy of the J'th column
                for(int i = 0; i < row; i++) {
                    try {
                        LUcolj.add(i, new Double(LU.get(i, j)));
                    } catch (MatrixIndexException ex) {
                        ex.printStackTrace();
                    }
                }
            
                for(int i = 0; i < row; i++) {
                    LUrowi = LU.matrix.get(i);
                    int kmax = Math.min(i, j);
                    double s = 0.0;
                    for(int k = 0; k < kmax; k++) {
                        s += LUrowi.get(k).doubleValue() * LUcolj.get(k).doubleValue();
                    }
                    LUcolj.set(i, new Double(LUcolj.get(i).doubleValue() - s));
                    LUrowi.set(j, LUcolj.get(i));
                }
            
                // Find pivot and exchange

                int p = j;
                for (int i = j+1; i < row; i++) {
                    if (Math.abs(LUcolj.get(i).doubleValue()) > Math.abs(LUcolj.get(p).doubleValue())) {
                        p = i;
                    }
                }
                if (p != j) {
                    for (int k = 0; k < col; k++) {
                        try {
                            double t = LU.get(p, k); 
                            LU.put(p, k, LU.get(j, k)); 
                            LU.put(j, k, t);
                        } catch(Exception e) {
                            e.printStackTrace();
                        }
                    }
                    int k = pivotArray[p]; 
                    pivotArray[p] = pivotArray[j]; 
                    pivotArray[j] = k;
                    pivsign = -pivsign;
                }

                // Compute multipliers.
                try {
                    if (j < row & LU.get(j, j) != 0.0) {
                        for (int i = j+1; i < row; i++) {
                            LU.put(i, j, LU.get(i, j) / LU.get(j, j));
                        }
                    }
                } catch(Exception e)  {
                    e.printStackTrace();
                }
            }
        }
    }
    
    public Matrix inverse() {
        this.LUDecomposition();
        System.out.println(LU);
        if(LU.row != LU.col) {
            return null;
        }
        return this.solve(new IdentityMatrix(this.row));
    }

    public Matrix solve(Matrix B) {
        // Copy right hand side with pivoting
        int nx = B.col;
        Matrix Xmat = B.getSubMatrix(pivotArray,0,nx-1);

        // Solve L*Y = B(piv,:)
        for (int k = 0; k < LU.col; k++) {
            for (int i = k+1; i < LU.col; i++) {
                for (int j = 0; j < nx; j++) {
                    try {
                        Xmat.put(i, j, Xmat.get(i, j) - Xmat.get(k, j) * LU.get(i, k));
                    } catch (MatrixIndexException ex) {
                        ex.printStackTrace();
                    }
                }
            }
        }
        
        // Solve U*X = Y;
        for (int k = LU.col-1; k >= 0; k--) {
            for (int j = 0; j < nx; j++) {
                try {
                    Xmat.put(k, j, Xmat.get(k, j) / LU.get(k, k));
                } catch (MatrixIndexException ex) {
                    ex.printStackTrace();
                }
            }
            for (int i = 0; i < k; i++) {
                for (int j = 0; j < nx; j++) {
                    try {
                        Xmat.put(i, j, Xmat.get(i, j) - Xmat.get(k, j) * LU.get(i, k));
                    } catch (MatrixIndexException ex) {
                        ex.printStackTrace();
                    }
                }
            }
        }
        return Xmat;        
    }
}
